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Abstract 

The initial distribution of gluons at the very early times after a high energy heavy 
ion collision is described by the bulk scale Qs of gluon saturation in the nuclear 
wavefunction. The subsequent evolution of the system towards kinetic equilibrium is 
described by a non-linear Landau equation for the single particle distributions [Q, |2|. 
In this paper, we solve this equation numerically for the idealized initial conditions 
proposed by Mueller, and study the evolution of the system to equilibrium. We dis- 
cuss the sensitivity of our results on the dynamical screening of collinear divergences. 
In a particular model of dynamical screening, the convergence to the hydrodynamic 
limit is seen to be rapid relative to hydrodynamic time scales. The equilibration 
time, the initial temperature, and the chemical potential are shown to have a strong 
functional dependence on the initial gluon saturation scale Qg- 
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1 Introduction 



An outstanding problem in high energy nuclear scattering is whether the hot and dense 
matter formed equilibrates to briefly form a plasma of deconfined quarks and gluons-the 
quark gluon plasma. This problem is of great topical interest, with collisions already taking 
place at the Relativistic Heavy Ion Collider (RHIC) and scheduled to take place several 
years hence at the Large Hadron Collider (LHC). 

Whether or not a quark gluon plasma is formed depends strongly on the highly non- 
equilibrium initial distributions of partons formed immediately after the collision. Clearly, 
these distributions must influence the subsequent evolution of the system towards equilib- 
rium. Furthermore, the problem is complicated by the rapid expansion of the system as 
a whole. The magnitude of the collision induced relaxation time relative to the expansion 
time is what determines whether equilibrium is indeed reached. 

The study of equilibration in relativistic heavy ion collisions is as old as the subject 
itself. However, very few ab initio studies exist that attempt to follow the evolution of the 
system, all the way, from the first instants of the collision to equilibrium 0, 1^, H . This was 
so because it was not known how to treat, in a self-consistent manner, the small Bjorken 
X "wee" parton modes that are responsible for particle production at central rapidities. 
These small x modes provide the initial conditions for the space-time evolution of the 
partonic matter formed in heavy ion collisions. 

Because the occupation number of small x modes in the nuclear wavefunctions is large, 
it was shown that classical methods could be used to compute their distributions 0. This 
classical effective theory is, in exact analogy to a spin glass, a color glass condensate 0, |[ , 
and is characterized by a bulk scale Qs - the momentum scale at which gluon distributions 
saturate. Analytical expressions for the parton distributions were obtained in Refs. |^ 
and |T^. Subsequently, quantum corrections to the classical "non-Abelian Weizsacker- 
Williams" fields were computed |TT[, and renormalization group methods ||^, |T^ devised 
to study how the classical parton distributions in the nuclei changed with energy (or 
equivalently, with x). The scale Qg now depends on x. This dependence is represented by 
a line in the x-Q^ plane - it separates the saturated, non-linear regime of QCD at high 
parton densities from that of linear QCD evolution. That Qs is a function of x, and grows 
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as one goes to smaller x, or equivalently, to higher energies, will only be implicit in this 
work. 

In the classical effective field theory approach, the problem of initial conditions can 
be formulated as the problem of finding solutions of the Yang-Mills equation with initial 
conditions given by the classical fields of each of the nuclei before the collision [|l^ . Since 
analytical expressions exist for the classical fields of the nuclei before the collision, the 
initial conditions are fully determined. 

Perturbative solutions of the Yang-Mills equations, which describe classical gluon pro- 



duction to lowest order, have been discussed by several authors |13, Ti . These were found 
to be infra-red divergent. Within the Yang-Mills approach, a fully non-perturbative treat- 
ment is therefore necessary. Non-perturbative, numerical solutions of the Yang-Mills equa- 
tions have been found recently In particular, one now knows the initial number 



and energy distribution of gluons after a collision |T^. An idealized form of the initial 
gluon distribution, in terms of the saturation scale of the color glass condensate Qs, was 
given recently by Mueller m. These idealized distributions are sufficient for the purposes 
of this paper. We will reserve a more quantitative analysis using the initial distributions 
of Ref. fl^ for a future work . 

The initial partonic system is completely out of equilibrium. The subsequent scattering 
and evolution of the system towards equilibrium was studied by Mueller [0, . In Ref. , 
he showed that, under the assumption that the system is undergoing boost invariant ex- 
pansion, the evolution of single particle distributions could be described by a non-linear 
Landau equation. At very early times, this equation can be linearized and studied an- 
alytically. However, the analytical approximations soon break down, and the evolution of 
the system cannot be followed analytically all the way to equilibrium. 

In this paper, we numerically solve the Landau equation proposed by Mueller. (Along 
the way, we compare our results to Mueller's analytical results and show quantitatively 
where the analytical approach breaks down.) We are thus able to follow the evolution 
of the system all the way to equilibrium. We study the dependence of the equilibration 
time, the initial temperature, and the chemical potential on the saturation scale Qs of the 
nuclear wavefunction, and on as- (Note: one may estimate ~ 1 GeV at RHIC and 
*to a Fokker-Planck equation 
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Qs ~ 2 — 3 GeV at LHC [|T], [T^.) We discuss the dynamical screening of the colhnear 
divergence arising from small angle scattering. Our results for equilibration are obtained, 



primarily, in a particular model of dynamical screening 20 



To conclude whether local kinetic equilibrium is attained in a system undergoing boost 
invariant 1+1-dimensional expansion, one should compare the equilibration time tequu to 
the hydrodynamic expansion time thydro- Typically, thydro ~ R/cq, where R is the radius 
of the nucleus, and cq is the speed of sound in the system. For an ideal ultrarelativistic 
gas, Co = 1/ -\/3. Our qualitative results suggest that for the energies of interest at RHIC 
and LHC, tequu "C thydro, thereby indicating that favorable conditions may exist for the 
formation of a hot gluon plasma at RHIC and at LHC. 

In studying the likelihood of equilibration, we have only considered number conserving 
2 — > 2 processes. Naively, 2 — > 3 processes will be suppressed by a power of as, and may 
be considered sub-leading. However, 2 — > 3 processes may, in principle, be much more 
efficient than 2 — > 2 processes in driving the system towards thermal equilibrium pT| . 
Whether this is indeed the case is a dynamical question we will not address in this work. 
Regardless, 2 — > 2 processes should set an upper bound on the equilibration time. Fur- 
thermore, it is likely that the numerical techniques developed here can be adapted to 
quantitatively study the effects of number changing processes. 

The paper is organized as follows. We begin in section 2 by discussing the initial 
conditions for gluon production in nuclear collisions. After the collision, particle production 
is described by the space-time evolution of classical gauge fields. At late times, when the 
system is dilute, one can define partons, and study their interactions in a kinetic approach. 
The kinetic transport equation derived by Mueller is described in section 3. In the following 
section, we discuss the problem of dynamical screening of collinear divergences, and describe 
the particular model of screening that is employed in this work. In section 5, we first 
briefly outline the numerical method, and compare results for a linearized (Fokker-Planck) 
version of the Landau equation with Mueller's analytical results. We next discuss numerical 
results from the solution of the full Landau equation. Various measures of equilibration are 
described and computed as a function of the saturation scale Qs and the coupling as- We 
summarize our results in the final section and comment on their experimental ramifications. 
We will also discuss the various uncertainties that need to be further quantified, and their 
likely impact on equilibration in high energy heavy ion collisions. An appendix contains a 
detailed description of the numerical method employed to solve the Landau equation. 
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2 The initial conditions for parton evolution 



The problem of initial conditions in heavy ion collisions is as follows. Before the collision, 
the nucleus, as viewed in a particular gauge, is a coherent superposition of various Fock 
modes containing differing numbers of partons |. Partons in Fock modes containing a large 
number of partons, \qqq ■ ■ ■ qgggg ■ ■ ■ qq), have small values of x-in the language of the par- 
ton model, the momentum of the nucleus is shared among a large number of constituents. 
Each parton has a + momentum k~^, and x = / , where is the momentum of 
the nucleus. The large multiplicities at central rapidities-which correspond to small x-are 
obtained when these small x Fock modes "go on-shell" . One should note that even though 
the small x modes carry small longitudinal momenta, their transverse momenta pt can be 
large. 

These small x partons are highly de-localized in the longitudinal x~ direction ^ - 
they have large wavelengths relative to the Lorentz contracted width of the nucleus. In 
addition, these virtual fluctuations are very short-lived. Any attempt to treat them in a 
naively classical transport approach is therefore problematic. An approach based on pQCD 
coUinear factorization is also inadequate unless the partons also have very large transverse 
momenta. This is because such a picture is predicated on convolving the probabilities of 
parton distributions in the two nuclei with the elementary parton-parton scattering cross- 
section. It cannot describe coherence effects which are important at small x. (For large pti 
these effects are suppressed.) 

It was realized some time ago that the small x Fock states in nuclei responsible for 

multi-particle production at central rapidities are states of high occupation number. As 

a consequence, their distributions in the nuclei are described by a classical effective field 

theory (EFT) The classical distributions in a single nucleus can be solved and an 

analytical form for the distributions obtained |]^, |TD[. The classical gluon distribution falls 

off as 1/ at large transverse momentum kt but saturates at smaller fc^-growing only slowly 

as \n{as^^/kt)/as- The infrared structure of the EFT is analogous to that of a spin glass 

in condensed matter systems [H-it is thus a colored glass condensate 0. 

tPor an illuminating exposition of this point of view, sec Ref. [ p^ . 
*if P+ is the momentum of the nucleus 
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The EFT approach was used in Ref. |T3] to treat the problem of initial conditions in 
nuclear collisions. It consists of solving the Yang-Mills equations in the forward light cone, 
after the collision, with initial conditions given by the (known) gauge fields of the nuclei 
before the collision. 

Gluons at central rapidities are produced classically. (Boost invariance is assumed.) 
Perturbative computations of the transverse momentum distributions of classically pro- 
duced gluons (lowest order in as and lowest non-trivial order in asQs/pt) were first com- 
puted in Ref. |[13[, and later by several authors [l^. In these computations, the number 
distribution is found to be infra-red divergent. They agree with the lowest order mini-jet 
predictions at large pt 0- 

When pt < Qs, all orders in Qs/pt contribute equally. The perturbative computations, 
to lowest order in Qs/pt are therefore not sufficient. Thus, even at lowest order in as, 
one needs to resum all orders in Qs/pt- An analytic solution to this non-perturbative 
problem has not been found. However, the problem was recently formulated as a 2+1- 
dimensional classical effective field theory and solved numerically (for an SU(2) gauge 
theory) on a 2-dimensional transverse lattice [0. The evolution of gauge fields is computed 



in real time, and at late times, the energy and the number |T3 of produced gluons at 
central rapidities can be computed. The initial distribution of gluons, from our numerical 
simulations is as follows. At low values of pt, pt < Qs it is remarkably like the Bose- 
Einstein distribution for a gas of massive particles in two dimensions. The gluons acquire 
a screening mass due to strong non-linear interactions - this also renders the distributions 
infrared finite. At large pt, Pt ^ Qs-, the distribution is a power law, Q'^/pj. The Born- 
pQCD prediction for the gluon distribution is therefore obtained, as expected, in the large 
Pt limit. 

The gluons produced at central rapidities are completely out of equilibrium. The study 
of the evolution of this system towards equilibrium was initiated by Mueller [|T], For 
simplicity, Mueller chose idealized initial conditions-a constant instead of a Bose-Einstein 
distribution for pt < Qg, and zero for pt ^ Q^-namely, the theta function 9{Q^ — pi). 
In this work, we will consider only these idealized distributions of Mueller. In a later 
work we will consider more realistic distributions computed very recently from the 
§For detailed comparisons, see the papers by Gyulassy and McLerran, and by Guo, in Ref. [n4[ 
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lattice simulations of Ref. |]T^. In general, the initial spatial distribution per unit rapidity 
is 

where is the transverse area, 77 is the space-time rapidity, and c, determined non- 
perturbatively, is a weak function of QsR- For an SU(2) gauge theory, it was computed 



recently in Ref. [[l^] to be c = 1.29 ± 0.09 in the regime of interest. 

The initial distributions are highly anisotropic, with the produced partons having zero 
longitudinal momentum p^- Assuming boost invariance, and using the relation 

(27r)3 dN 

for the single particle distribution f{x,p), one finds at late times that [] 

fix,p) = -^jS{p.)e{Q',-pl). (3) 
asNc t 

At late times, the occupation number of partons in the transverse plane becomes small, 
and their evolution can no longer be followed on the lattice. However, when the occupation 
number is not too large, transport theory can be applied to study the further evolution of 
the system. The initial condition for this evolution is given by the single particle distribu- 
tion in Eq. ^. 

Using the results of Ref. [ffH] , which determined the time at which gluons come on-shell 



after a collision, we can estimate the time at which small angle scattering between the 
gluons can be described by the Boltzmann equation. The formation times, using Qs = 1 
GeV for RHIC and Qs = 2-3 GeV for LHC, are 

U ~ 1.40 GeV"^ for RHIC , U ~ 0.62 GeV^^ for LHC (4) 

Next, we suppose that a transport theory based treatment of gluon scattering is applicable 
when the gluon occupation number is less than unity. Letting drj he dt] = dz/ti = 1 or 

dN 



dz = ti, and using rf^a, from 0] at 2; = 0, we find 



^0 = (5) 

Here to > ti is the time where the gluon occupation number is dilute enough that their 
subsequent interactions can be described as the small angle scattering of on-shell gluons. 

^The (5-fimction distribution is of order l/Ap^- We shall see later, from the early time analytical 
solution of the Fokker-Planck equation, that the longitudinal momentum changes little in the time in 
takes for the occupation number to decrease below unity. 
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3 The Landau transport equation 



At times t > Iq, where to was defined in Eq. ^, a Boltzmann-like transport equation is 
appropriate for describing tlie late-time evolution of the highly anisotropic initial gluon 
distribution discussed in the previous section. We will show below that, for central high 
energy heavy ion collisions, assuming boost invariance, the transport equation reduces to 
a Landau-type transport equation . 

In the problem of interest, the typical scale of spatial variations is large compared to 
the typical scale of gluon-gluon scattering. The evolution of the system can therefore be 
described by a local Boltzmann-like kinetic equation for the single particle gluon distribu- 
tion f{x,p, t), which gives the density of excitations of momentum p at a point (x, t). The 
Boltzmann equation is 

^ + v.Vf = C{f), (6) 

where C{f) is the collision integral which represents the change, due to collisions, per unit 
time, in the number of particles per unit phase space volume. 

Let us first consider the left hand side of the above equation. We will assume that the 
transverse dimensions of the collision volume in central heavy ion collisions are sufficiently 
large that the initial expansion of the system is one-dimensional. The distribution / 
then only depends on z, the co-ordinate corresponding to the collision axis. Since it is 
independent of the transverse co-ordinate Xt, Eq. ^ has the form 

^^^vJ-^^ = Cm,z,t)], (7) 

where Vp^ = Pz/\p\- A further simplification follows from the assumption that the central 



rapidity region is approximately Lorentz invariant under boosts As was pointed out 
by Baym many years ago [0], this assumption greatly simplifies the problem of solving 
Eq. since it relates the distribution function at different z's in the central region. It 
therefore suffices to compute / in the zero-rapidity slice alone. Because / is a scalar under 
Lorentz transformations, it satisfies the relation 

f{PuPz,z,t) = f{pt,p',,t) . (8) 

Here p'^ = 'j{pz — up), the transformation velocity u = z/t, and 7 = (1 — n)^^/^ = t/r, 
where the proper time r = 



z 



2 
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Computing df /dz using Eq. the Lorentz transformation relation, dT/dz\z=Q = 0, 
and dp'z/dz\z=o = —p/t, one finds at 2 = 0, 

""^dz ~ t dp,- ^ ' 

The Boltzmann equation now reduces to 

Note that, as a consequence of our manipulations, assuming uniform distributions in the 
transverse direction and boost invariance in the longitudinal direction, the single particle 
distribution f{x,p, t) is now expressed as a function of p and t alone. The kinetic equation 
above, in the absence of collisions, has the free streaming solution for / |^^, namely. 



f{t,p) = exp (-/5oVPi+pSV7f) 



where (3q and 70 are constants. 

We shall now consider the collision term C[f] on the right hand side of the above 
equation. If we assume that the changes in momentum g in a collision are small, namely, 
q/Qs ^ 1, we can treat the collision integral as simple diffusion in momentum space. 



Following the discussion in Lifshitz and Pitaevskii [iq], one can write 



C[f] = -P^, (12) 

OPa 



where the flux 



Note that v = pl\p\- 

We have used in Eq. |13| the cross-section for gluon-gluon elastic scattering, which in 
the limit of small momentum transfer squared = — t, is 

da _ fasN.^ i2iTy 1 

di \ 71 J 2(iV2 - 1) £2 • ^ ) 

The collision integral contains a logarithmic coUinear divergence arising from small angle 
scattering. It is represented here by C, deflned as 

C= ^. (15) 



The minimal momentum transfer is, gmin = QsGmm, where 6mm is the minimum scattering 
angle. Also, gmax is the maximum momentum transfer. How one may determine q^i^ and 
gmax is discussed in the next section. 

Having started from Eq. ^, combining Eq. and Eq. |T^, we arrive at the considerably 
simpler expression, the Landau equation |jl8|, ^ 

df Pz df ^ ^ 2 



dt t dpz 

with the definitions 



XnC Vif + 2A n_i C ■ (vf), (16) 



^y't)=9G j7^f{t,p) , n^iit) = gc , \ = 2nas—^ -. (17) 



Here, go = '^{N^ — 1) and Vp- refers to differentiation with respect to p. 
Multiplying both sides of Eq. ^ by d^p and integrating, one finds 

|m=0, (18) 

namely, the number density of gluons has the behavior n oc 1/t. Taking the second moment 
of Eq. |1^, one finds the following exact relation in the central region |^ 

where the local energy density is given by 

e{t)=gG J -^,PfM^ (20) 
and the longitudinal pressure is 

P.W=9c/|^f/(p,*). (21) 

In the hydrodynamic regime where Pl = e/3, the ideal gas equation of state gives e ~ t~^^^. 
We will return to a discussion of these quantities later on in the paper. 

Following Mueller, it is convenient to define the quantities r] = tn and ?7_i = tn_i. The 
Landau equation can thus be re-written as 

^dn^_^^dfMl. = Ar/£ V|/(t,p-) + 2Ar/_i(t)£ V,^- {vf{t,p)). (22) 

If ?7_i is taken to be independent of time, this equation is the Fokker-Planck equation: 
the coefficient of the first term would be the diffusion coefficient, while that of the second 
would be the coefficient of friction. 
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4 Screening of collinear divergences in the Landau 
transport equation 



In the one gluon exchange approximation, as in Coulomb scattering, the transport integral 
in Eq. |l^has a logarithmic collinear divergence-represented by £. In Ref. [Q, the condition 
applied at early times was that in the maximal distance 1 / gmin corresponding to the minimal 
momentum transfer gmiru one has at most one scattering f\ From this condition, by using 
the initial gluon distribution, Mueller finds at the earliest times that 

A more general expression, for any arbitrary number density , is easily derived to be 



Alternately, one can use the more conventional prescription of regulating the logarithmic 
collinear divergence in C through the exchange of dressed gluons characterized by a Debye 
screening mass. A kinetic theory expression for the Debye screening mass can be derived 



from the screening of a time-like gluon propagator in a medium of gluon excitations |^ 
One obtains 

^ f d^pM__^^. V,V(t,p). (25) 

vr^ \q\^oJ q-p 

One can check that for an equilibrium Bose distribution, one recovers the standard result: 
= AnasNcT"^ /3. For our initial conditions, it is the screening of the longitudinal 
gluons in the transverse direction that is relevant. Performing an integration by parts, this 
transverse Debye mass is given by the expression 



At the very early stage of the 2^2 scattering, the expression in Eq. ^ is, in principle, 
less reliable since the system is completely out of equilibrium. However, at later times, as 



II This condition satisfies the requirement that the soft field (due to aU the other hard gluons) that a 
particular hard gluon scatters off be at most of size 1/g. This is so because the freed gluons come from 
the saturated component of the nuclear wavefunction of size 1/g-see the discussion in section 2. 
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a = 0.3, c = 1 .3, Q. = 1.0 (GeV) 




Figure 1: The dynamical screening masses, Mm from eq. (^) obtained from ||T1| and the 
mx from eq. ( PB| ) obtained from hnear response theory plotted as functions of time in 
GeV~^ units. Also plotted versus time are the logarithms Cm. and Ct obtained by combining 
Eq.'s (0) and respectively with Eq. (^. 

the system approaches equilibrium, m£, should be more reliable. One could attempt to 
parameterize the infrared cut-off of small angle scattering such that it interpolates between 
the two limits. However, at present there is no theoretical justification of any particular 
form. Since collisions become more frequent as one approaches equilibrium, it is likely that 
the Debye mass is a more reliable cut-off to use. At any rate, it is important that one at 
least have the right limit as the system approaches equilibrium. 

In Refs. [|l], 0, the value of gmax is held fixed at early times. This is because if one 
re-writes C in Eq. [15| as 

£ = In f ^ V (27) 
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it is argued, on the basis of analytical approximations to Eq. |^ below, that ^ changes 

1 

appreciably only over a long time scale {t ~ e^^). Since the Debye mass has a stronger 
time dependence, C grows with time. On account of the approximations made in arriving 
at this conclusion, it is useful to check whether this behavior of C is obtained by computing 
(pt) and rriT self-consistently at each step in our numerical simulations 0. We indeed do 
find, as shown in in Fig. |], that C grows rapidly initially as a function of time but levels 
out to a constant as the system approaches equilibrium . 

In Fig. |l|, we plot the infrared masses in Eqs. and as well as the corresponding 
contribution C to the collision integral, as a function of time (in units of 1/Qs). One 
sees that while the relative difference at early times is large, it is much less so at later 
times. We have also solved the transport equation using both cut-offs. While there is some 
quantitative difference between the two, they do not affect our qualitative conclusions. For 
the rest of this work, we will therefore use the Debye mass in Eq. ^ 

Debye screening occurs due to the exchange of longitudinal gluons. One may worry 
whether the screening of transverse gluons will introduce an additional scale since the static 
magnetic mass is parametrically larger than the Debye mass. However, the conventional 
wisdom is that, for transport cross-sections, magnetic screening is also cut off in the infrared 



by the Debye mass [^, The same, for instance, is not true of color transport-the color 



conductivity is regulated by a magnetic mass of order g^T The dynamical screening 



of infrared divergences in 2 — n (n > 2) transport cross-sections merits further study. 

5 Results from numerical solution of the Fokker- 
Planck Equation 

We will first consider, as a "warmup exercise", the solution of Eq. ^ for the case where 
?7_i is held constant. The Landau equation then reduces to the Fokker-Planck equation. 
To check the consistency of our numerical simulation, we will compare, in detail, the 
numerical solution of the Fokker-Planck equation to the early time analytical solution 
found by Mueller in Refs. ^ and p|. We will discuss solutions to the non-linear Landau 
equations in the next section. 



*We thank A. Dumitru for discussions on this point. 
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Following Mueller, we first define / = tf, C = |log(t/to) = where to is a constant 
defined in Eq. ^. (As discussed in the previous section, this assumption is justified at early 
times.) Scaling A/3 A, we can rewrite ( p2l) as 



" Wz^') ^^^"^^ + ''-^ ■ ^^^^ • ^^^^ 

Eq. (^) is a Fokker-Planck equation if ?7_i is held fixed. If one integrates Eq. ^ with 
respect to d^p, d^pp\ and d^ppl we obtain Eqs. (17), (18) and (19) in Ref. 

d , , 

+ 2{p^) = 2A,?-^ftK), (30) 

and 

|(pi> = 4A,4l-!tl|±l)+^P,K). (31) 



where Pl = tPi is defined by eq. ([2l|). These equations are exact. We can therefore 
perform the following consistency check on our solution. With the initial conditions in 
Eq. H and the boundary condition that f{px,Py,Pz,t) = if any p^, Py or pz oo, we 
numerically determine / for all time, directly compute and (p^), and check to see if 
Eqs. (|29|) to (0) are satisfied. Mueller supposes that at early times the first term on the 
l.h.s and the last term on the r.h.s of Eq. ^ are small, thereby yielding the solution 

{pD = Xvt (32) 

It is demonstrated in Fig. |^ that this approximation, in practice, is not a very good 
one. At very early times, if ,^ ~ we see that the first term on the l.h.s of Eq. ^ is of 
the same order as the remaining terms. We cannot therefore ignore it. We can ignore the 
term proportional to Piit) in Eq. ( |5DD for early times. Then Eq. ( |5D| ) becomes: 

d 



^^{pD + 2{pI) = 2Xr^C 



Its solution is 
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Figure 2: Numerical solution to the Boltzmann equation with constant rj^i where each 
term in Eq. ( P^ is plotted as a function of ^ = ln(t/to). 

{pl) = Xv^ + ^{exp{-20-l). (33) 

This solution (see Fig. ^ agrees much better with the numerical results than Eq. ( p2D 
does. We define Cj{^) as 

2Ar/e - ^Pl{0 - UpD - ^{pD = C/(0, (34) 
where C/(^) should be small if Eq. (|30[) is satisfied. Fig. ^ verifies indeed that this is 

so. 

Next, look at Eq. (0). If we assume that both ?7_i and r^+i are constant with respect 
to time and determined by the initial condition Eq. (^), 
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Using Eq. (|35| ) we re-write Eq. (|3TD as: 

ii(,i)._*,,«H-^p,(0. (36) 

Again the term containing Pl is small. We therefore solve the equation above, with the 
appropriate boundary conditions, to find 

(pi) = f - l^r^e. (37) 
As previously, define Cii{C,) such that 

^(pD + - ^pLiO = CniO- (38) 
at, 3 r] 

Fig. 1^ verifies that Eq. ( pH]) is satisfied for early times. 

Finally, we check to see how well Eq. ( pQ]) is satisfied numerically. Fig. ^ plots the 
numerical result for t] vs. t and demonstrates that rj is constant with respect to t. As 
discussed earlier. Fig. |^ also compares the numerical result for {pD and {p]_) with the 
analytical solutions of Eqs. (^) and ( PP] ) given by Eqs. ( P^D and 



6 Results from numerical solution of the Landau 
transport equation 

In the previous section, we numerically solved the Boltzmann equation using Mueller's 
approximation for the minimum scattering angle, and with ?7_i fixed for all time (Eq. 
(p8|)) Since the previous analysis is only expected to be valid for early times, one must 
solve the general expression ( |22D with ?7„i determined self consistently. 

In this section, we will discuss results from numerical solutions of the Landau equation 
in Eq. ^ employing the initial single particle distribution given by Eq. ^ and the dynamical 
screening mass given by Eq. The equation is a second order partial integro-differential 
equation. It can be solved by combining the Alternating Direction Implicit (ADI) method 
'I''!' Recall that is an integral representing the moment of /. 
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Figure 3: Each term in Eq. ( ^8|) plotted as a function of ^ = \n(t/to] 



and the Crank-Nicholson differencing scheme |3^ . A detailed discussion of the numer- 
ical procedure can be found in Appendix A. 

We will first begin by briefly describing the evolution of the single particle distribution 
and its moments. Next, we will discuss various measures of kinetic equilibration and 
hydrodynamic flow for a system undergoing boost-invariant one dimensional expansion. 
Finally, we will discuss our numerical results for these quantities for different values of Qs 
and asf^- 



6.1 Single particle distributions 

The initial distribution for the numerical solution of Eq. (p2| ) is given by (^, at the initial 

time, determined from Ref. [Q, given by Eq. (^). The time evolution of the single particle 

ttWe will ignore running coupling effects in this analysis. It is argued in Ref. that running coupling 
effects in the evolution are 0(^^05) and therefore suppressed. 
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Averages, Nc = 3, c = 2.5, a = 0.2 




0.3 - 
0.2 - 
0.1 - 




time 

Figure 4: Numerical results for (p^) and plotted as a function of time, compared 

with the analytical approximations in Eqs. p3| ) and ( |57D . Also plotted is 77 as a function 
of time. 

distribution f{t,Pz,p±) is shown in Fig. ^for a particular set of initial parameters. It begins 
its time evolution as a delta function in (represented practically by a narrow Gaussian 
distribution) and as a step function in p± — see Eq. (^). As time proceeds, we see the p± 
modes get scattered off the transverse plane in the longitudinal direction. The larger p± 
modes decrease rapidly as the pz distribution widens out. At about t = 11.95 GeV~^, the 
Pz distribution is at its broadest extent. 

The behavior of the single particle distributions are seen more clearly in Fig. where 
we have plotted /, as a function of Pz or p^, for fixed p^ and pz respectively. We see initially 
that the two distributions are widely different. At time proceeds, they converge relatively 
rapidly at soft momentum but much more slowly in the tails. If thermal equilibrium is 
defined strictly as the distributions being completely isotropic, this condition is reached 
only asymptotically, if at all. 

Fig. 1^ shows, for particular initial parameters, the averages (p^), {pI) and {pD as 
functions of time. {pD starts at zero and quickly rises before converging slowly to (p^.) and 
(Py), which, in turn, decrease monotonically with time. The non-trivial behavior of {pD 
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f(t,p,,p2) for different times (GeV ^) 




Figure 5: The number of gluons per unit of phase space f{t,p) plotted as a function of 
p± and Pz for different times. Nc = 1, as = 0.3, c = 1.3 and Qg = 1 GeV. Units of time 
are in GeV~^. 

is because the system is undergoing longitudinal expansion. In a box at rest, one expects 
that {pD will show more of a monotonic behavior before leveling off. It again appears from 
Fig. 1^ that the convergence to isotropic distributions is very slow. This is particularly so 
since the second moment of the distribution weights the high momentum tail unduly. The 
latter, as we observed in Fig. |, takes longer to equilibrate. 

Interestingly, as we will discuss in the following, the convergence of bulk thermodynamic 
observables to the expected equilibrium values is much more rapid. A likely explanation is 
that these observables are much less sensitive to the high momentum tail of the distribution. 
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Single particle distributions 




1.5 



p, or p. (GeV) p, or p. (GeV) 



Figure 6: The single particle distributions, for fixed px or pz plotted as a function of pz 
and Px respectively, for different times in the evolution of the distribution. The results are 
for as = 0.3, Qs = I, c = 1.3 and Nc = 3. 

6.2 Kinetic equilibrium and hydrodynamics 

From Figs. and (|^), it appears that the distribution, as a whole, become isotropic 
only asymptotically. This statement is particularly true of the tails, in pz and pt, of the 
distributions. The distributions agree more closely at softer momenta. Nevertheless, as 
we shall discuss below, when we look at thermodynamic signatures of equilibrium, they 
converge relatively rapidly to the expected behavior in a fluid undergoing boost-invariant 
1-dimensional expansion. 

Let us first discuss what this expected behavior is. Since we only consider 2^2 
processes in this work, the total number of gluons is fixed. The equilibrium solution of the 
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Figure 7: The averages (p^.), (p^) and (p^) versus time in GeV ^ units for different values 
of asi c and Qg. In all cases, = 3. 

classical Boltzmann equation, Eq. (0), thus has the form 

f{t,p) = exp[f]{t){fi{t)-p)]. (39) 

Here fi{t) is the chemical potential, and /3{t) = 1/T(t), where T is is the temperature. 
Substituting this equilibrium distribution in Eq. (pIT]), one obtains 

n(t)=2M^r^(t)exp(|) , n_i(t) = ^^^2^nt)exp(^) . (40) 

The average energy density is defined to be 

e{t)=9aj-0^,\p\f{t,p). (41) 

From Eq. (^) and Eq. (^), the energy per particle at equilibrium is, 

E{t) = e{t)/n{t) = 3T{t). (42) 
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Now, the entropy density of a classical Boltzmann gas is defined as 



^W = -^c/(^/log/. (43) 
In equilibrium, the entropy per particle is simply, 

5(t) = 4l = 3-^. (44) 

We noted previously, since the number of gluons is conserved, that tn{t) = constant. 
Also note that, since the entropy per particle is constant in equilibrium, we find from the 
above equation that fi/T = constant. From these constraints, the system in equilibrium 
must satisfy 

t = constant , e t'^^^ = constant . (45) 
Finally, recall we had defined the longitudinal pressure in the central slice in Eq. ^ as 



(27r)3 p 

One can similarly define the transverse pressure Pt to be 

Sp p1 



From Eq. 19, and the above relations, the condition for ideal hydrodynamics is 



PT = PL = \e. (47) 

The approach to equilibrium, in the sense of "saturating" the above thermodynamic (and 
hydrodynamic) identities, has been studied, in the relaxation time approximation, by sev- 



eral authors 32, R3]. This approximation is also employed in studies with mini-jet 



initial conditions ||31|, The collision kernel in Eq. [1^ can be written formally as [p5| , |26[] 

as 

C[f{p..Put)] = -^-tlp^ ^ (48) 

where 6', the collision or relaxation time, is in general a function of time and momentum. 
The relaxation time approximation is one where the momentum dependence of 9 is ne- 
glected. Baym studied equilibration in this approximation taking ^ to be a constant. He 
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showed that the thermodynamic relations in Eqs. H3] and were satisfied only asymp 



totically in time Q Nevertheless, the convergence of the system to the asymptotic value is 
rapid. The system, while not quite in local thermodynamic equilibrium, is sufficiently close 
to it that equilibrium is a good working assumption. Subsequently, Gavin studied equili- 
bration in the relaxation time approximation assuming 9 to have the time dependence |^ ; 



9 = at, where a is a constant. Depending on the value of a, the system approaches the 
hydrodynamic limit quickly (smaller a's) or free streaming (larger a's). Heiselberg and 
Wang [Q, studied the general case, 9 = t^. They conclude that for p < 1 thermalization 
is attained, while the system free streams for p > 1. The case p = 1 studied by Gavin is 
the marginal one, interpolating between the two regimes for different a's. In a subsequent 
paper Heiselberg and Wang study the dependence of 9 in ffnite temperature QCD, 



and tentatively conclude that p ~ 1/3 < 1. 

6.3 Convergence to the hydrodynamic Umit 

In this work, as discussed in section 3, we go beyond the relaxation time approximation in 
treating the collision kernel. It is not necessary therefore, once a screening mechanism is 
postulated, to specify the time dependence of the collision time. In the following, we will 
discuss our results for the thermodynamic relations stated in the previous section. 

In Fig. ^ we show the entropy per particle, S, plotted as a function of time. As shown for 
typical values of as, c and Qs, it initially increases rapidly and ffattens out, monotonically 
approaching its asymptotic value (thereby implying that fi/T goes to a constant as well). 



The expression in Eq. used in the computation of Fig. |g, is the correct one, except at 

very early times when the well-known full expression for the entropy density of a Bose gas 

should be used. Thus, where the result in Fig. § is going to zero is where the classical 

Boltzmann expression for the entropy density is breaking down. We have checked that 

the full expression ensures that the entropy per particle is always positive definite. The 

triangles in Fig. ^ correspond to the entropy per particle assuming the partons are initially 

localized on the two dimensional transverse plane. 

The convergence of other thermodynamic quantities to the hydrodynamic limit is shown 

in Fig. ^ Again, we note that the convergence to their asymptotic values is much more 

*This is indeed what we would conclude from our study of the single particle distributions in the 
previous sub-section. 
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rapid than one would expect by looking at the single particle distributions alone. Our 
results suggest that the collision time effectively has a time dependence with p < 1. 




Figure 8: The entropy per particle S plotted as a function of time for some typical values 
of as, c and Qs- In both cases, Nc = 3. The triangles denote the entropy per particle of a 
two dimensional Boltzmann gas. 

Since the convergence to the hydrodynamic limit is only asymptotic, deciding when the 
system can be described in terms of thermodynamic quantities is somewhat subjective. 
Here, we define the equilibrium time teg as the time it takes for Tt^^^ (|i2| ) and the entropy 
per particle S (|4^) , to reach 90% of their maximum asymptotic value (see Figs. ^ and |). 
One can then also extract the "initial" temperature and chemical potential that correspond 



to teg by using Eqs. ^ - ^ 



The equilibration time tgg is very sensitive to the values of a^, c and Qg. Table | shows 
the equilibration time t^g, the initial temperature T^g, and the chemical potential fi for 
typical values of as, c and Qs- The parametric behavior of these quantities (for a fixed 
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value of the non-perturbative constant c) as a function of Qs-, for two different values of 
the coupling constant as-, is also shown in Fig. 0. Larger values of as and Qs yield smaller 



teq- 



a = 0.3, c = 1 ,3, Q. = 1 .0, Nc=3 
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Figure 9: The top panel plots Tt^^^ = t^^'^ where e is the energy density and n the 
number density; the middle panel the e, the longitudinal pressure Pl and the transverse 
pressure Pt] the bottom panel e^^^ and the entropy density s; all plotted versus time in 
GeV~^ units, for typical values of as, c, Qs and N^.. 

These results can be qualitatively understood as follows. One can show (see Eq. 
that the equilibration time is parametrically 

~ k (vS) ■ 

In Fig. 1^, we note that teq decreases roughly as 1/Qs- Also, it is greater for smaller as 
as one would expect from this expression. Similarly, from requiring that tn{t) = constant. 
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one finds for the initial temperature tliat 



We see tliat tliis dependence on Qg and as is confirmed in Fig. p!o| . 
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Figure 10: The estimated values of the time, teg, temperature Teq, and the chemical po- 
tential, /i, at the onset of equilibrium, plotted as functions of Qs for as = 0.2 and 0.3. 

At RHIC, we expect (roughly) that as ~ 0.3, c ~ 1.3 and Qs ~ 1.0 GeV. For the 
idealized initial conditions discussed here, the corresponding time and temperature at ki- 
netic equilibrium are teg ~ 3.2 fm and T^g ~ 174 MeV. At LHC, we expect ~ 2-3 
GeV, therefore teg ~ 2.4 fm - 1.4 fm and Teg ~ 321 MeV - 472 MeV. Are these numbers 
realistic? Not so (especially at RHIC), since it is unlikely that, at these temperatures, 
the system is a weakly coupled gluonic gas. At the time scales and temperatures corre- 
sponding to a rapid convergence of the system to the hydrodynamic limit, other (likely 
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as 


Q.(GeV) 




re,(MeV) 


/i(MeV) 


0.1 


1.0 


12.24 ± 1.92 


171.05 ±4.9 


137.78 ± 19.3 


0.2 


1.0 


5.8 ± 0.83 


166.28 ± 5.19 


144.71 ± 15.22 


0.2 


2.0 


4.48 ± 0.75 


304.32 ± 19.39 


212.95 ± 31.23 


0.3 


1.0 


3.24 ± 0.44 


174.27 ± 5.21 


157.86 ± 16.18 


0.3 


1.4 


2.77 ± 0.60 


234.70 ± 11.60 


195.71 ±31.75 


0.3 


2.0 


2.36 ± 0.44 


320.72 ±22.38 


249.61 ±25.06 


0.3 


2.4 


1.80 ± 0.56 


386.38 ±28.48 


337.56 ±26.18 


0.3 


3.0 


1.42 ± 0.45 


471.69 ±36.49 


457.72 ±36.17 



Table I: Equilibration times, temperatures and chemical potentials as a function of as and 
Qs for c = 1.3 and = 3. 

non-perturbative) effects might become important. One cannot conclude definitively that 
this is the case because we have not considered realistic initial distributions nor have we 
discussed the importance of particle number changing processes. We will comment on this 
point in the final section below. 

7 Summary and outlook 

We have solved numerically a non-linear transport equation, Eq. (|2^) , which describes 
the evolution, after a heavy ion collision, of single particle gluon distributions in the cen- 
tral rapidity slice. The initial conditions for the solution of this equation are the highly 
anisotropic, idealized, initial conditions first discussed in the context of parton transport 
by Mueller [Q, |^. The distributions are controlled by a single scale-the saturation scale 
Qs of parton distributions in the nuclear wavefunction before the collision. Only elastic 
gluon-gluon scatterings are treated. Equilibration, in this approach, is dominated by small 
angle scattering. The coUinear divergences that occur are regulated dynamically by a cut- 
off of the Debye form. We have checked that, at early times, the linearized Fokker-Planck 
equation reproduces the analytical results of Mueller. These analytical approximations 
however cannot be carried through to times relevant for equilibration. 

We find that the tails of the initially highly anisotropic initial conditions converge very 
slowly to the expected isotropic equilibrium distribution. This behavior is confirmed by the 
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behavior of {pD and (p^). Despite the slow convergence of the single particle distributions 
to the isotropic thermal shape, thermodynamic observables such as the entropy per particle, 
the energy per particle, and the transverse and longitudinal pressures converge more rapidly 
to the hydrodynamic behavior expected of 1-dimensional, boost invariant expansion. The 
more rapid convergence for these observables occurs because they are more sensitive to 
softer momentum modes and less so to the high momentum tails. 

Using a particular criterion for equilibration (that thermodynamic observables have 
reached 90% of their asymptotic value), we extracted the equilibration times and the initial 
temperatures and chemical potentials. We have studied how they behave for varying values 
of Qs and as- Even though small angle scattering is very inefficient (the equilibration time 
is an order of magnitude greater than the formation time), it is nevertheless smaller than 
the hydrodynamic time scale thydro ~ -R/cq- Here R is the radius of the nucleus and cq 
is the speed of sound in the fluid. The relatively long equilibration times correspond to 
relatively low temperatures, T <^ Qg. At these temperatures, it is unlikely that the system 
can be described as a weakly coupled gluon gas. It is therefore reasonable to ask whether 
other effects, be they perturbative or non-perturbative in nature, may significantly alter 
our results []. We will enumerate below, in order of increasing complexity, those effects that 
are amenable to a weak coupling treatment. 

Firstly, recall that our results were obtained for idealized initial conditions. More 



realistic initial distributions have recently become available [|1^], and they are qualitatively 
different from the idealized distributions. Secondly, since the occupation numbers of the 
gluons are large initially, the effects of Bose enhancements in the Boltzmann equation 
should be taken into account. Both of these effects are straightforward to incorporate in 
our approach, and results with these improvements will be reported shortly [p^ . 

Thirdly, as we discussed in the introduction, number changing 2-^3 processes may be 
very important in establishing kinetic as well as chemical equilibrium |]SB|. Even though 
these processes are suppressed by an additional power of as, they are more efficient in 
re-distributing the momenta. Unfortunately, past treatments have been handicapped by 
uncertainties in how one treats infrared divergences. Work in progress suggests that this 
problem may be cured, and that equilibration is much more rapid ||2l|. If this is indeed 
tPor an interesting recent take on this topic see Ref. pS. 



28 



the case, then initial temperatures are relatively close to the initial saturation scale, and 
the weak coupling treatment of equilibration will have been self-consistent. 

Ultimately, the goal is relate the rich variety of hadronic and electromagnetic spectra 
that will soon be available at RHIC (and some years later, at LHC) to properties of the 
initial nuclear wavefunction, that may independently be probed in deeply inelastic eA or 
in pA collisions. This paper is a first quantitative step in that direction. 
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A Numerical analysis of second order partial integro- 
difFerential equations 

Eq. (^) is a second order partial integro-differential equation whose numerical solution 
requires some care to ensure stability. Fortunately, standard finite differencing schemes 
can be used and numerical stability guaranteed. Eq. (p2D, is of the form: 

where u = u(t, x, y, z). 

To discretize Eq. (|49|) we first define a mult i- dimensional grid: 

t = to + n^T , X = xo + «A , y = yo + j A , z = zq + kA, 

where n = 0, 1, 2, ... A^, i = 0, 1, 2, ... J, j = 0, 1, 2, ... J, and A; = 0, 1, 2, ... /sT. Next we 
employ the Alternating Direction Implicit (ADI) Method using the Crank-Nicholson finite 
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differencing scheme The ADI method is especially useful in solving parabolic multi- 
dimensional equations on rectangular grids. For problems with three spatial dimensions, 
such as Eq. (^91), the ADI method is implemented by splitting each time step of size AT 
into three steps of size AT/3. At the current fractional time step, n + ^ for example, only 
one of the spatial derivatives are evaluated and the others are evaluated at the previous 
time step, n. If we choose not to split the time step in such a way, the solution to Eq. (^), 
after discretizing, would require us to invert a large matrix of the same size as our grid. 
The purpose of splitting the time step in such a way is to reduce the problem in question 



into the solution of a tridiagonal matrix (as shown below). To illustrate this, Eq. ( ^Qf ) 
is discretized using the ADI method with the Crank-Nicholson finite differencing scheme. 
First discretize and evaluate the x derivatives at the current time step n + ^: 



_ n 
AT/3 



A2 



"^i+ljk ~^ '^i-ljk 2u.;i,^^ + U 



1 
'3 
-Ijk 



n 

ijk "jj+lfc "ij-lfc 
~'^'Kjk + ^ij-fc+l + '^Ijk-l 

+ linear terms. 



(50) 



n + 



Increment the time step by 1/3 and evaluate the y derivatives at current time step 

2. 



'^ijk" 



u 



n+i 



AT/3 



A2 



u. 



i+ljk 



-Ijk 



'^"'ijk 



2 
3 

ij+lk 



■"'3 "-^3 
~'^^ijk + ^ijk+1 + 



U. 



ij—lk 



■^^ijk 



-\- linear terms. 



(51) 



Increment the time step again by 1/3 and evaluate the z derivatives at the current time 
step n -|- 1: 



u. 



n+l 
ijk 



u 



ijk 



AT/3 



A2 



n+l n+l 
'^i+ljk + '^i-ljk 



'^^ijk 



I n-t-g 

+ '^ij+lk 



I "+3 



2 
3 

ij—lk 



.2i/"+3 , „,n+l , n+l 



-|- linear terms. 



(52) 
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In Eq.'s (^Uj) to (0), "linear terms" refers to the terms in (^) which contain first order 
derivatives. Since the Crank-Nicholson finite differencing scheme is stable only if the terms 
with the highest order of derivatives are implicit, the lower order terms are allowed to be 
either explicit or implicit. We assume here for brevity in notation that the linear terms 
are explicit and evaluated at the previous time step. Therefore, we do not write the linear 
terms out. 

The solution to Eq.'s ( pOD to (|52[) at every time step, is simply the solution of a tridi- 
agonal matrix at every n. For example, Eq. (|50D can be arranged as follows: 



«"«ri5. + + T-^r+ifc = ^.--fc, (53) 

where 
and 



en 



AnAT 



ijk 



+ linear terms. 



Eq. ( p3|) is simply an equation of the form: 

M ■ u = d, (54) 

where M is a tridiagonal matrix. Since we know M and d at every previous time step 
n, we can in principal solve for u at every current time step + |, given the boundary 
conditions u° Ugj^, u"ofc and u^q. The inversion of sparse matrices (such as tridiagonal 
matrices) is usually numerically trivial. 

In solving Eq. (^) one needs to specify the boundary conditions suitable to solving 
a second order differential equation. In this work we have specified u^jf. by Eq. (|^) and 
required that u vanish at the boundary of space: u^j^, u^q^, u^jq, u]j^, u'^jf, and m^j^ = 0. 
Furthermore we require that the first derivatives vanish at x,y,z ±oo. Therefore, 
Mjjfc — ~ 0, Uiji^ — Mojfc ~ 0, and so on for j and k. 
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